/* Release 3.10 */
/*-------------------------------------------------------------------------*
*                                                                          *
*  COMPLEX.C - complex number class implementation.       		   *
*                                                                          *
*  COPYRIGHT (C) 1989..1992 Clarion Software Corporation.                  *
*  All Rights Reserved                                                     *
*                                                                          *
*--------------------------------------------------------------------------*/

#include <clib.h>
#include <complex.hpp>

#pragma module(init_prio=>22)

double norm(const complex& c) {

  double r, i;

  r=real(c);
  i=imag(c);
  return r*r + i*i;
}

double arg(const complex& c) {

  return atan2(imag(c), real(c));
}

complex operator*(const complex& c1, const complex& c2) {

  double ires, rres;
  double r1, r2, i1, i2;

  r1=real(c1);
  i1=imag(c1);
  r2=real(c2);
  i2=imag(c2);
  rres = r1*r2 - i1*i2;
  ires = r1*i2 + r2*i1;
  return complex(rres, ires);
}

complex operator/(const complex& c1, const complex& c2) {

  double ires, rres;
  double r1, r2, i1, i2;

  r1=real(c1);
  i1=imag(c1);
  r2=real(c2);
  i2=imag(c2);

  rres=(i1*i2+r1*r2)/(r2*r2+i2*i2);
  ires=(r2*i1-r1*i2)/(r2*r2+i2*i2);

  return complex(rres, ires);
}

const complex& complex::operator*=(const complex& c2) {

  double ires, rres;
  double r2, i2;

  r2=real(c2);
  i2=imag(c2);
  rres = re*r2 - im*i2;
  ires = re*i2 + r2*im;
  re=rres;
  im=ires;
  return *this;
}

const complex& complex::operator/=(const complex& c2) {

  double ires, rres;
  double r2, i2;

  r2=real(c2);
  i2=imag(c2);

  ires=(im*r2 - i2*re)/(i2*i2 + r2*r2);
  rres=(im-r2*ires)/i2;
  re=rres;
  im=ires;
  return *this;
}


double abs(const complex& c) {

  double i, r;

  i=imag(c);
  r=real(c);
  return sqrt(r*r + i*i);
}

complex acos(const complex& c) {

  complex r;

  r=sqrt(1 - c*c);
  r=complex(-imag(r), real(r) );
  r=log(c + r);
  return complex(imag(r), -real(r));
}

complex asin(const complex& c) {

  complex r;

  r=sqrt(1 - c*c);
  r=log(complex(-imag(r), real(r)) + r);
  return complex(imag(r), -real(r));
}

complex atan(const complex& c) {

  complex r, t;

  t=complex(-imag(c), real(c));
  r=log((1.0+t)/(1.0-t))/2.0;
  return complex(imag(r), -real(r));
}

complex cos(const complex& c) {

  complex t;

  t=complex(-imag(c), real(c));
  return (exp(t) + exp(-t))/2.0;
}


complex cosh(const complex& c) {

  return (exp(c) + exp(-c))/2.0;
}

complex exp(const complex& c) {

  double ex;

  ex=exp(real(c));
  return complex(ex*cos(imag(c)), ex*sin(imag(c)));
}

complex log(const complex& c) {

  complex t;

  t=arg(c);
  t=complex(-imag(t), real(t));
  return log(abs(c)) + t;
}

complex log10(const complex& c) {

  complex t;

  t=arg(c);
  t=complex(-imag(t), real(t));
  return (log(abs(c)) + t)/log(10.0);
}

complex pow(const complex& base, double expo) {

  return exp(expo*log(base));
}

complex pow(double base, const complex& expo) {

  return exp(expo*log(base));
}

complex pow(const complex& base, const complex& expo) {

  return exp(expo*log(base));
}

complex sin(const complex& c) {

  complex t;
  complex twoi(0.0, 2.0);

  t=complex(-imag(c), real(c));
  return complex((exp(t) - exp(-t))/twoi);
}

complex sinh(const complex& c) {

  return (exp(c) - exp(-c))/2.0;
}

complex sqrt(const complex& c) {

  complex t;

  t=sin(arg(c)/2.0);
  t=complex(-imag(t), real(t));
  return sqrt(abs(c))*(cos(arg(c)/2.0) + t);
}


complex tan(const complex& c) {

  return sin(c)/cos(c);
}

complex tanh(const complex& c) {

  return sinh(c)/cosh(c);
}

ostream& operator<<(ostream& st, const complex& c) {

  st<< '(' << real(c) << ',' << imag(c) << ')';
  return st;
}

istream& operator>>(istream& st, complex& c) {

  double r=0.0, i=0.0;

  if(st.peek() == '(')
    st.get();
  st >> r;
  if(st.peek() == ',')
    st.get();
  st >> i;
  if(st.peek() == ')')
    st.get();
  c=complex(r, i);
  return st;
}